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Abstract 

Diffusion phenomena in a multiple component lattice Boltzmann Equation 
(LBE) model are discussed in detail. The mass fluxes associated with differ- 
ent mechanical driving forces are obtained using a Chapman-Enskog analysis. 
This model is found to have correct diffusion behavior and the multiple diffu- 
sion coefficients are obtained analytically. The analytical results are further 
confirmed by numerical simulations in a few solvable limiting cases. The LBE 
model is established as a useful computational tool for the simulation of mass 

transfer in fluid systems with external forces. 
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I. INTRODUCTION 



The lattice Boltzmann Equation (LBE) method is an increasingly popular method of 
computational fluid dynamics. As an extension of the lattice gas cellular automaton [|l],3, 
this method simulates fluid motion by following the evolution of a prescribed Boltzmann 
equation instead of solving the Navier-Stoke equations. At the macroscopic level, it has been 
proved that the Navier-Stokes equations can be recovered from the Boltzmann equation. 
There have been many publications on this subject, and interested readers are referred to 
these publications |^,^ and the references therein for the history, background and details of 
this method. Recently, convincing numerical simulations have shown that the LBE method 
can simulate fluid flow at high Reynolds number with very good accuracy PJ^. 

An important advantage of the LBE method is that, since it deals with the distribution 
functions, microscopic physical interactions of the constituent fluid particles can be conve- 
niently incorporated. For complex fluid flows with interfaces between multiple phases and 
phase transitions, the complex macroscopic behavior is the consequence of the interactions 
between the fluid particles. Since the early stage of the development of the lattice gas and 
lattice Boltzmann method, considerable effort has been invested in incorporating particle 
interactions into the lattice models so that complex fluid behavior including multiphase 
flows can be simulated. Rothman and Keller |^ developed the flrst lattice gas model for two 
immiscible fluids. A Boltzmann version was formulated later In this scheme, the particle 
distributions of the two species are rearranged in the interfacial region in a way dependent 
on concentration gradients. The same idea was also used to reduce the diffusivity in a mis- 
cible two-component system [0]. Flekk0y introduced another two-component LBE model of 



two miscible components |[11|J10|| , in which the sum of the distribution functions of the two 
components and the difference of them are made to relax at difference rates to the speci- 
fled distribution functions so that the diffusivity is independent of the viscosity of the fluid 
mixture. In another lattice gas model of liquid- vapor phase transition [Q, the long-range 
interaction was introduced by exchanging momentum over several lattice spacings. 
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In a previous publication we presented an LBE model for multiple component 



systems which includes interactions between particles of the same and different components. 
An interaction potential is defined for each of the components, and an additional momentum 
exchange is introduced as the consequence. By considering nearest-neighbor interactions 
only, we were able to alter the equation of state of the fiuid to a general class of functional 
form, allowing the simulation of non-ideal gases and their mixtures. With this model, we 
can simulate the motion of the interfaces and mass transfer between different phases. The 
components in the system can be completely miscible or partially immiscible depending 
on the temperature and the relative strengths of the interactions. Given the interaction 
potentials, we have also analytically obtained the coexistence curve, the density profile 



across a liquid- vapor interface, and the surface tension . 

In many real-world multiphase problems, mass transfer in the presence of external forces 
is involved. An example is the centrifugal separation of components of a fiuid mixture. In 
centrifugal separation, a large acceleration is applied to the fiuid mixture and the particles 
of different components diffuse through the mixture at different rates. Moreover, when 
dealing with problems involving phase transitions such as dissolution and evaporation, the 
boundaries of the flow fleld usually move. Because the lattice Boltzmann method treats 
the multiphase problem using a single equation, many complicated effects can be naturally 
integrated into the algorithm. 

Before the LBE model can be used to quantitatively simulate complex fluid flow problems, 
we must have a thorough understanding of the behavior of the model itself. It is essential 
that the physical parameters of real systems can be matched. In a previous publication 



T5| , we derived the macroscopic equations for the multiple component lattice Boltzmann 
model and obtained the mutual diffusivity in a binary mixture by calculating the decay rate 
of an inflnitesimal concentration perturbation. We found that the diffusivity depends on 
the collision times, the concentrations of the components, and the interaction potentials in 
a complicated way. The diffusivity can be tuned to be arbitrarily close to zero and even 
negative. 



In this paper, we provide a detailed study of diffusion in the multiple component LBE 
model including interparticle interactions and external forces. Compared with other LBE 
miscible model, the current scheme has the following features. First of all, the diffusion in 
this model is Galilean invariant. The diffusivity is independent of the flow velocity. Sec- 
ond, multiple diffusion in a system with arbitrary number of components can be simulated. 
Third, since the constituent components can be either miscible or immiscible to each other 
depending on the interaction force, this model can be used to simulate diffusion in a multi- 
phase system with mass transfer between different phases. In Section |l|, we briefly review 
the multiple component non-ideal gas lattice Boltzmann model and give the macroscopic 
fluid equations satisfied. In Section |ITTt mass fluxes are calculated using the same Chapman- 
Enskog technique that we previously developed. The effects of different mechanisms which 
drive the diffusion are identified and the multiple diffusion coefficients are obtained. The 
results are further confirmed by numerical simulation of a few analytically solvable limiting 
cases. In Section |V|, we give the conclusion and offer some more discussion on simulation 
of multiple component systems with this model. 



II. THE MULTI-COMPONENT LBE MODEL 

For completeness, we briefly review the LBE model for a multiple component fluid mix- 
ture with interparticle forces ||13H15|. Consider a lattice gas system in D-dimensional space 
with particles of S components moving from one lattice site to its b nearest neighbors and 
colliding with each other at lattice sites at each time step. The particles of the ath com- 
ponent have the molecular mass, m„. The distribution function of the particles of the crth 
component is assumed to evolve according to the following Boltzmann equation: 

<(x + e„, t + 1) - <(x, t) = -- k(x, t) - nl^^^^ (x, t)] (1) 

where, {Sq; a = 1, . . . , 6} is the set of vectors of length c pointing from x to its h neighbors; 
n^(x, t) is the population of the particles of component a having velocity at lattice site x 
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and time t. The collision term on the right-hand side takes the form of single-relaxation-time 
for each component. This collision term has the BGK form, named after Bhatnagar Gross 
and Krook |T^. It can be efficiently implemented on computers. It has been shown that 
with a properly chosen distribution function, ?t,J^^'^'^^(x, t), the correct Navier-Stokes equation 
can be recovered from Eqs. (|I|) at the macroscopic level [0,|r^,0. 

For the multiple component LBE model, we chose n^^^'^^ = Naijic^, u^'^), where ria = J2a 
is the number density of the crth component. The functional form 

Na{n,n) = <^ ^ (2) 
[n[d^-^y, a = 

is the same one that yields the correct Navier-Stokes equations for a single component 
LBE model. For simplicity, we chose = do in previous publications. This choice is not 
required to obtain the Navier-Stokes equations. We will allow d^ to be different for for each 
component. The parameters, u^*^, in the above distribution function are chosen to be 

p^u'^ = p^u' + r^F,, (3) 

where p^ = m^na is the density of the ath component and Fg- is the total external force 
acting on particles of the ath component. F^. includes both external forces and interparticle 
forces. With u^^ so chosen, at every site and for every collision step, each component gains 
an additional momentum F^ due to external and interparticle forces. In the absence of any 
additional forces, all the components are assumed to have a common averaged velocity u'. It 
follows from the requirement that the total momentum must be conserved at each collision 
when Fg- = that 

u' = Ef^E<e.)/X:^. (4) 

(T = l \ 'o" a=l J (T = l ^o" 

In general, this averaged velocity is different, and should be carefully distinguished, from 
the fluid velocity which represents the overall mass transfer. 

In previous publications p^JT^ , we incorporated an interparticle force between the par- 



ticles at sites x and x'. This interparticle force is proportional to the product of a function 
of the particle number densities. 



s 

F.(x) = -V^.(x) ^ Yl ^-(^' x')^^(x')(x' - x) + p^g., (5) 

x' 5-=l 

where, ^^^^(x, x') is the Green's function which satisfies Q^a = Qua, and ipai'^) = V^o-(^(x)) 
is a function of the number density which plays the role of an interaction potential. The 
form of this function directly determines the equation of state, as will be seen later, g^. is 
the external force acting on the ath component, and it can be different for each component. 
If only nearest-neighbor interactions are included, the above expression becomes 

s b c% ^ 

Fa = -i^a + ^<^)^<^ + P'^S'^ - — JT^^ Q^^^i^^ + P^Sa- (6) 

a=l a=l ^ a=l 

With the forces included, the sum of the momenta of all the components is not 
conserved at each site by the collision operator, although the total momentum of the whole 
system is still rigorously conserved. Since the macroscopic fluid velocity of the mixture, 
u, represents the over-all mass transfer rate, it should be deflned by the total momentum 



averaged before and after each collision |T5|. A Chapman- Enskog expansion procedure can 
be carried out to obtain the following momentum equation for the fluid mixture as a single 
fluid: 

— -f u ■ Vu = ^ + Y ^-g- + '^V'u, 7 

^ P ^1 

where, p = Z]f=i Po- is the total density of the fluid mixture, and = pa/p is the mass 

fraction of component a. The pressure p is given by the following non-ideal gas equation of 

state: 



c2 



^(1 - d„)m^n^ + ^ XI ^fraV'aV'a 



(8) 



Here, the flrst term on the right-hand-side is the kinetic contribution and the second term 
is the potential contribution due to interparticle interaction. Notice that for a mixture of 
ideal gases, the partial pressure does not depend on the molecular mass of each component. 
To make the equation of state @ approach that of a mixture of ideal gases in the limit of 
very weak interactions, it is appropriate to chose 1 — = Dcl/rjiaC^, where cq is the sound 
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speed in the mixture in the absence of interactions. Eq. (|^) in this case can then be written 

as 



(9) 



Without losing generahty, we treat the da^ as arbitrary constants. 

With each component having a distinct d^r, and with the additional terms due to the 
external force paga, the mass conservation equation that we derived before for component 
a has to be modified slightly: 



dpa_ 
dt 



+ V ■ (p^u) 



1 



c\l-d^) 
D 



(T=l 



(10) 



It can be verified that, when summed over all the components, the right-hand-side of 
Eq. ([T0|), which represents the diffusion of the component a, is zero, so that the conti- 
nuity equation for the whole fluid is satisfied. Generally, in the presence of large-scale fluid 
motion, the diffusion of the components through each other is coupled to the large-scale 
flow. The evolution of each component is governed by the most general macroscopic fluid 
equations (0) and (|^). However, in many cases the fluid is at rest except for the motion 
caused by the diffusion of the different components. We will discuss in detail the diffusion 
in a fluid mixture in next section. 



III. DIFFUSION IN THE MULTI-COMPONENT LBE MODEL 

The components of a fluid mixture are said to be diffusing into each other if the mean 
velocities of the components differ. The local velocity of the fluid mixture can be defined 
in several different ways: by averaging the velocities of the constituents by mass, mole 
or volume [|1^]. The diffusion velocities are then defined relative to this local velocity. 
Mathematically, all the averaging methods are equally useful in describing the diffusion 



of the constituents. Following the treatment in Chapman and Cowling pOf, we use the 



mass fluxes of the components in our calculation. Here again, since the momentum of each 
component changes at each collision, to obtain the over-all mass flux, we must average before 
and after collisions (cf. ref. ll^). We have 



E 



(11) 



By applying the same Chapman-Enskog technique previously developed [0, at the second 
order, namely n"^ = n'^'^^^ + n'^^^\ we obtain the relative mass flux of the crth component 
after tedious but straightforward manipulations: 



D 



cr=l 



:i2) 



The velocity u^- — u is the mass-averaged diffusion velocity of component a with respect to u, 
which indicates the motion of component a relative to the local motion of the fluid mixture. 
Noticing that the right-hand-side terms of Eqs. (|1^) are exactly what the divergence operator 
acts on in the right-hand-side of Eq. (|10|), we can simply re- write Eq. (|T0|) as the continuity 
equation of the crth component 



dt 



+ V ■ {Pa^a) 



0. 



(13) 



We have therefore demonstrated that each component satisfies its own continuity equation 
at the second order. 

Following the convention in the diffusion literature, we define the mass flux of component 
cr as jo- = Pcr(uo- — u). Obviously we have X^o- jo- = 0. From Eqs. ([l^), we can attribute the 
generation of the mass flux j^- to three different driving mechanisms: the concentration 
gradients, the pressure gradient, and the inequality of the external forces acting on different 
components. The diffusions driven by these driving mechanisms are called the ordinary 
diffusion, pressure diffusion, and the forced diffusion respectively. It is convenient to separate 
the effects of the different diffusions from each other by decomposing the mass flux into its 
corresponding parts 
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i _ ]{x) I "(p) I '(3) 



(14) 



where j^'^^ j^^-*, and j^^^ are respectively the mass fluxes associated with ordinary, pressure 
and forced diffusions. The mass flux corresponding to the forced diffusion, j^^-*, can be 
separated most easily from the others by collecting in Eqs. (|^) the terms containing go-. 
The result is 



-per ^{Xi - Si„)Ti igi-Yl ^kgk 
i=l \ k=l / 



(15) 



where 5ia is the usual Kronecker delta. The mass flux of component a depends on forces on 
all the components. It can be shown that the sum of j^^^ over all the components vanishes. 
If all the gj are the same as in the case where gravity is the only external force, all the j^^^ 
vanish. Therefore, forced diffusion only occurs when the external forces applied to all the 
components are not equal. An example of such a system is a mixture of differently charged 
particles in external electric field. In addition, when all the are equal, the j^^^ can also be 
simplified: 



After the separation of j^^^ , the terms remaining in Eq. 
combination of the density gradients with the help of Eqs. ( 



(16) 

can be written as a linear 
and (I), 

(17) 



where the coefficients D^ji are 



+ E 

fc=i 



k=l 



Qkii'k'4'i- 



(18) 



Here, the first term is the ideal-gas contribution. The second term is the potential part 
due to interactions. Since variations in both the mass fraction and the pressure will cause 
density variations, to separate the effects of the mass fraction variation and the pressure 
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variation, we must write the mass flux in terms of the gradients of the mass fraction and 
the pressure. Using the deflnition pi = pxi, Eq. (|1^) can be written as 

i'a^ + j^"^ = ^ (pEi^^^Vx, + Vp|:z^..x,^ . (19) 

By taking the gradient of Eq. (|]), we have 

c% ^ c% ( ^ ^ \ 

= 7T E A- Vp, = 77 P E Vx,- + Vp E (20) 
^ i=i ^ \ i=i i=i / 

where Aj = (1 — dj)/h + ijj'^ Yl,i=i Qiji^i- EUminating Vp from the two equations above, we 
obtain 

j=i \ l^j=i^jXj I 2^j=i^jXj 

The mass fluxes associated with concentration gradients and the pressure gradient can be 
immediately identifled as 



■,2u^ S 



Jcr p. 

^ 1=1 



S IS 

E(-Dai^j — D^jAi)xj I E ^j^j 



Vx,, (22) 



^^^=Vpj:D^,xJj2A,x,. (23) 

i=i / j=i 

The coefficients in front of the mass fraction gradients are the multiple diffusion coefficients 
of our LBE model. Noting that ^iXi = 1, the mass flux of component a can be written as 
being dependent on the mass fraction gradients of all but the ath component. 

Using Eq. (pTSP and (p^)-([23|), we computed the contributions to the mass flux from 
the three mechanical driving forces. Except for thermal diffusion, this LBE model has the 



correct types of diffusion behavior compared with the continuum theory of diffusion |]T9 



The reason for the lack of thermal diffusion is that this current model assumes that the 
temperature is a constant and independent of space. 

Ordinary diffusion given by Eq. (|22|) has rather complicated dependence on the gradients 
of all the concentrations. We have analytically given the multiple diffusion coefficients in 
terms of the interaction potential ipi, the collision interval Xj, the constants 1 — di which 
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defines the mole volume of each components in ideal gas limit, and the mass fractions Xj. 
The multiple diffusion coefficients are concentration-dependent, and can be theoretically 
adjusted to simulate specific material properties. 

This LBE model also exhibits a pressure diffusion phenomenon. Depending upon the 
parameters, when a pressure gradient is applied to the mixture, there could be net fluxes 
of individual components in an originally homogeneous mixture. We can therefore use this 
model to study problems such as centrifuge separation. When different external forces are 
applied to the individual components in the mixture, an originally homogeneous mixture 
will separate so that concentration gradients will be generated to balance the effect of the 
forced diffusion. While all these different types of diffusion occur in the LBE system, the 
mixture satisfies the Navier-Stokes equations as a single fluid. 

The diffusion mass fluxes given by Eqs. ([15|) and p^)-(p3D are for the most general case 
and are valid for systems with an arbitrary number of components and for arbitrary forms of 
interaction potentials, as long as the interaction is not so strong that the components become 
immiscible and segregate into different phases. Since the coefficients have a complicated 
dependence on the parameters and the densities themselves, analytical solutions of the 
densities are generally difficult to obtain even in the "static" case. We will discuss a few 
limiting cases in which the density distributions can be analytically solved and compare the 
results with those from numerical simulations. 

A. Diffusivity in a binary mixture 

We will derive the mutual diffusivity in a binary mixture (5* = 2) using the mass flux 
obtained above and compare it with previous results. For a binary fluid mixture, Fick's first 
law of diffusion can be stated in our notation: 

jl = -pVWxr, (24) 

which gives the definition of the mutual diffusivity T>. Since xi + X2 = 1, and Vxi = — Va;2, 
after some tedious algebra, Eq. (|22|) can be written in the form of Fick's law, with the 
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diffusivity, 



If we set dcr = do, this can be easily verified to be identical to the mutual diffusivity we 
obtained previously |T5| by computing the decay rate of an infinitesimal concentration per- 
turbation. This result has been confirmed by measurement of the actual decay rate of a 
concentration wave in numerical simulations with the LBE code [|l5l . The derivation here is 
more general in the sense that no linearization of equations is required. 



B. Mixture of ideal gases 

A mixture of ideal gases can be simulated by setting Q^-a = 0, and choosing the constants 
di so that 1 — di = Dc^/niiC^. In this case the pressure is proportional to the total number 
density of the mixture. The second term in Eq. (p!8D vanishes, and = (1 — di)/b. If the d^ 
are all equal, we can verify using Eq. (|T8]) that J2i D^iXi = 0, and therefore the mass fluxes 
j^P) vanish identically. This implies that for an ideal gas mixture, pressure diffusion occurs 
if and only if the components of the mixture have different molecular weights. 

We consider the case in which a common conservative external force, given by g = — V0, 
is applied to all the components. Forced diffusion does not occur in this situation and the 
condition for equilibrium is j^^'-* + j^^-* = 0. By directly substituting into Eq. ([l^), we can 
confirm that the following density profiles satisfy the equilibrium condition: 

(26) 

where are constants determined by the initial conditions. 

Figure ^ shows the steady-state density profiles in a two-component numerical simulation 
performed on a two-dimensional hexagonal lattice with 16 sites in the x-direction and 256 
sites in the ?/-direction. Due to the effect of the non-square lattice, the actual dimension 
in lattice units is 16 x 128-\/3. A periodic boundary condition is used in the x-direction. 



Per = P„ exp 



c\l-d^) 
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In the ?/-direction, a solid wall is placed at x = 0, and bounce-back boundary conditions 
are used at the wall. The constants, di and 6/2, are chosen to be 0.4 and 0.6 respectively. 
The external force potential is chosen to be = ~g{y/LY without losing generality, where 
g = 0.1 is a constant and L is the y dimension of the lattice. Plotted are the typical measured 
density profiles at equilibrium, together with their theoretical solutions given by Eq. 
The agreement is always excellent independent of the parameters such as Xj and the mean 
densities. 



C. Forced diffusion 

The effects of forced diffusion is also examined for a binary mixture of ideal gases. The 
constants are chosen to be equal for the two components {di = do) to eliminate the effects 
of pressure diffusion. The mass fluxes of the two components in this case consist of the 
contributions of ordinary diffusion and forced diffusion. The steady-state density profiles of 
the two components are now given by the equations j^^-* + j^^-* = 0. This equation can be 
simplified if we assume Ti = T2 = r. We have 

- VVXi + XiX2T{gi - g2) = 0, (27) 

where V is the mutual diffusivity, in this case, ^-^^^^(r — 1/2). gj = — V</'i is the external 
force acting on component i. Clearly when gi = g2, no forced diffusion can occur and the 
steady-state mass fraction profiles are homogeneous. With gi 7^ g2, we can solve Eq. 
to obtain 



Xi 

— = Ci exp 

X2 



(28) 



V 

where Ci is an integration constant determined by the over-all mass ratio of the two com- 
ponents. In Figure 0, we confirmed this solution by numerical simulation with the same 
geometry and boundary conditions as before, except that now we have two different force 
potentials which are chosen to be 0i = gi{y/L) and 02 = 92{y/LY, with gi = —0.1 and 
92 = 0.1. 
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IV. CONCLUSION 



In this paper, we discussed in detail the diffusion behavior in a previously proposed 
multiple component LBE model. The effects of particle interaction and external forces are 
included in the analysis. We calculated, using the Chapman- Enskog expansion, the mass 
fluxes in the mixture due to different driving mechanisms and we obtained the multiple 
diffusion coefficients. The LBE model is found to exhibit all types of diffusion except the 
thermal diffusion. All types of diffusion are Galilean invariant. The analytic calculation is 
consistent with numerical simulations in several solvable limiting cases. 

With the diffusion coefficients analytically calculated and the effects of external forces 
identified, we are now able to quantitatively simulate a wide class of practical problems 
involving diffusion, separation and fluid flow simultaneously. After the transport phenomena 
are satisfactorily treated, chemical reactions among components can also be added easily in 
this model to simulate chemical reaction processes. 

Since diffusion in a multi-component fluid is itself a very complicated phenomenon, the 
calculation of the transport coefficients from the parameters of the model is tedious but 
straightforward. For practical engineering applications, this process can be automated. 

Finally, we point out again that since this model only simulates isothermal fluids, the 
possibility of directly simulating an interesting thermal diffusion phenomenon, known as the 
Soret and Dufour effects is ruled out. We consider this as an important area for improvement 
in this LBE model. 
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FIGURES 

FIG. 1. Equilibrium density profiles of the two components in a binary mixture of ideal gases 
with a pressure gradient appUed. The theoretically predicted profiles are plotted as solid lines and 
the numerical results are plotted as symbols. 

FIG. 2. Density ratio of the two component in a binary mixture of ideal gases at equilibrium. 
Two different external forces are applied to the two components separately. The solid line is the 
analytical solution and the symbols are numerical results. 
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